Efficient simulation of relativistic fermions via vertex models 
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We have developed an efficient simulation algorithm for strongly interacting relativistic fermions 
in two-dimensional field theories based on a formulation as a loop gas. The loop models describ- 
ing the dynamics of the fermions can be mapped to statistical vertex models and our proposal is 
in fact an efficient simulation algorithm for generic vertex models in arbitrary dimensions. The 
algorithm essentially eliminates critical slowing down by sampling two-point correlation functions 
and it allows simulations directly in the massless limit. Moreover, it generates loop configurations 
with fluctuating topological boundary conditions enabling to simulate fermions with arbitrary peri- 
odic or anti-periodic boundary conditions. As illustrative examples, the algorithm is applied to the 
Gross-Neveu model and to the Schwinger model in the strong coupling limit. 

PACS numbers: 02.70.-c, 05.50.+q, ll.10.Kk, 11.15.Ha 



Simulating strongly interacting fermions, like in Quan- 
tum Chromodynamics (QCD) or in Nambu-Jona-Lasinio 
models, is considered to be rather difficult and continues 
to be a challenge. Contrary to common lore, this is not 
directly due to the fact that the fermionic degrees of free- 
dom are Grassmannian variables, but rather due to the 
non-locality of the determinant which is obtained upon 
integrating out the fermionic fields. Moreover, simula- 
tions of fermions arc usually hampered by critical slow- 
ing down towards the chiral limit where the fermions be- 
come massless and the correlation length of the fermionic 
two-point function diverges. The established standard 
method to perform such calculations on the lattice is to 
use the Hybrid Monte-Carlo algorithm [l[ which deals 
with the non-locality of the determinant by rewriting 
it as an integral over bosonic " pseudo-fcrmion" fields. 
The algorithm then requires to deal with the inverse of 
the fcrmion Dirac operator, however, the operator be- 
comes ill-conditioned towards the massless limit and the 
simulations slow down dramatically. In this letter we 
propose a novel approach which circumvents the above 
mentioned problems. It is based on a (high-temperature) 
expansion of the fermion actions which reformulates the 
fermionic systems as g-state vertex models, i.e., statisti- 
cal closed loop models. In particular, the method is di- 
rectly applicable to the Gross-Neveu (GN) model and to 
the Schwinger model in the strong coupling limit. These 
models can be shown to be equivalent to specific vertex 
models @, H, 0, Q and our simulation method, based on 
a proposal by Prokof'ev and Svistunov is effectively a 
very efficient updating algorithm for generic vertex mod- 
els (in arbitrary dimensions). In fact, the algorithm es- 
sentially eliminates critical slowing down and is able to 
simulate the fermionic systems at the critical point and 
directly in the massless limit. 

We start with illustrating the reformulation in terms of 
closed loops in the GN model. The model is most natu- 
rally formulated by employing Major ana fermions @, @|. 
Here we are using Wilson's Euclidean lattice discretisa- 



tion for which the action density of the model is 
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where £ is a real, two component Grassmann field de- 
scribing a Majorana fermion with mass to, C = —C T is 
the charge conjugation matrix, and d,d*,d denote the 
forward, backward and symmetric lattice derivative, re- 
spectively. The Wilson term ^d* d, responsible for re- 
moving the fermion doublers, explicitly breaks the dis- 
crete chiral symmetry £ — > 75^, £ T C — > £, T C-f5 and re- 
quires a fine tuning of to — > m c towards the continuum 
limit in order to restore the symmetry. A pair £1 , £2 
of Majorana fermions may be considered as one Dirac 
fermion using the identification ip = ^(£1 + 4> = 
-^(£f — i£j)C and the corresponding GN model with 
TV Dirac fermions has an 0(2N) flavour symmetry. At 
g = 0, integrating out the Grassmann variables yields the 
partition function in terms of the Pfaffian 



^GN = Pf 



C(7/A + m- 2 d * d ) 



2N 



(2) 



For g 7^ one usually performs a Hubbard-Stratonovich 
transformation and introduces a scalar field a conjugate 
to £ T C£ which yields the new fermionic action 

X 

-X^^ci^Or + A), (3) 

together with an additional Gaussian Boltzmann factor 
exp{— l/(2g 2 ) J2 X &{x) 2 } for the scalar field. 

In order to reformulate the model in terms of closed 
loops (or equivalently dimers and monomers) we follow 
the recent derivation of Wolff (see 0,(1 for alternative, 
but more complicated derivations). One simply expands 
the Boltzmann factor for the fermionic fields and makes 
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use of the nil-potency of the Grassmann variables upon 
integration. Introducing tp(x) = 2 + m + cr(x) and the 
projectors P(±/x) = (l=F7u)/2 we can write the fermionic 
part of the GN path integral (up to an overall sign) as 



i(x) 



n(c T wcp(^(x+A)) 



(4) 
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FIG. 1: The vertex configurations and weights of the eight- 
vertex model (top row) and the extended model (bottom row). 



where m(x) = 0, 1 and bp,(x) = 0, 1 are the monomer 
and bond (or dimcr) occupation numbers, respectively. 
Integration over the fermion fields yields the constraint 
that at each site m(x) + \ ^ M b^{x) = 1. Here the sum 
runs over positive and negative directions and b-^(x) = 
b^(x — fi). The constraint ensures that only closed and 
non-intersecting loops of occupied bonds contribute to 
the partition function and also accounts for the fact 
that the loops are non-backtracking, a consequence of 
the orthogonal projectors P(±^i). The weight uj(£) of 
each loop I can be calculated analytically [§| and yields 
= 2~ nc -/ 2 where n c is the number of corners along 
the loop. The sign of ui will generically depend on the 
geometrical shape of the loop [9( prohibiting a straight- 
forward probabilistic interpretation of the loop weights 
in dimensions d > 2. 

In two dimensions, however, the sign of the loop only 
depends on the topology of the loop, as recently clar- 
ified by Wolff 0, and is determined by the fermionic 
boundary conditions. It is therefore useful to classify 
all loop configurations into the four equivalence classes 
£oO: £io, An, £11 where the index denotes the total wind- 
ing (modulo two) of the loops in the two directions. The 
weights of all configurations in £10 and C\i for exam- 
ple will pick up an overall minus sign if we change the 
fermionic boundary condition in the first direction from 
periodic to antipcriodic, while the weights of the config- 
urations in Cqq and £01 remain unaffected. As a con- 
sequence, if we sum over all the topological equivalence 
classes with positive weights, i.e., Z = Zc 00 + Zc 10 + 
Zc 01 + Zc ±1 we effectively describe a system with un- 
specified fermionic boundary conditions. Vice versa, the 
partition function Z|° = Z Caa + Z Cia - Z Cai + Z Cll , e.g., 
describes a system with fermionic b.c. antipcriodic in the 
first and periodic in the second direction. 

Before describing our method to simulate the loop for- 
mulation of the GN model, essentially generating closed 
loop configurations according to their loop weight, it is 
useful to point out the equivalence to the 8- vertex model 
13, El- The model s formulated in terms of the eight 
vertex configurations shown in the top row of Fig. [T]with 
weights u>i,i — I,... 8. The partition function is then 
defined as the sum over all possible tilings of the square 
lattice with the eight vertices such that only closed (but 



possibly intersecting) paths occur. To be precise, one has 



CP x 



(5) 



where the sum is over all closed path configurations (CP) 
and the weight of each configuration is given by the prod- 
uct of all vertex weights in the configuration. For the GN 
model we have the following weights 



u> 2 



0, 



U> 3 

UJ 5 = UJ 6 = UJ 7 



U>4 



1, 



^ = 72' 
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i.e. each corner contributes a factor l/v2, while crossings 
of two lines are forbidden (uj2 = 0) and each empty site 
carries the monomer weight u)\ = tp(x). From here it also 
becomes clear that for a single Majorana fermion the in- 
teraction term proportional to g is irrelevant. Since the 
partition function is now factorised into terms at each x 
we can integrate over the scalar field a(x) at each site sep- 
arately. However, since the integration measure is even 
in <j(x) the term linear in cr(x) will not survive and the 
monomer weight reduces to lj± = 2 + m. This is how the 
free Majorana fermion is recovered after the Hubbard- 
Stratonovich transformation. If one considers two or 
more Majorana flavours which are coupled through the 
four-fermion interaction, the integration over u(x) be- 
comes non-trivial, but can still be done analytically Q 
yielding 8- vertex models coupled to each other. It should 
be stressed, however, that from an algorithmic point of 
view the generic case with an arbitrarily varying field 
cr(x) is equally accessible and involves no complication 
whatsoever. 

The fact that the GN model with a single Majorana 
fermion is effectively a free fermion system expresses it- 
self also through the vertex weights fulfilling the free 
fermion condition 

.10 

vertex models fulfilling the free fermion condition as well 
as those in zero field 3, 14 1 are analytically solvable, in 
particular also the standard Ising model. We use this fact 
to our convenience and use the Majorana GN model as 
a benchmark which allows to compare the results of our 
algorithm with analytic results. Instances of the 8-vertex 
model for which no analytic solutions are known, but are 
accessible with our algorithm, include the Ising model 
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with additional next-to-ncarcst-ncighbour and quartic in- 
teractions . Another 8-vertex model with a fermionic 
interpretation is the one-flavour Schwinger model with 
Wilson fermions in the strong coupling limit @, EH . The 
vertex weights are given by 



LJl = (to + 2) 2 , 

UJ 2 = 0, 



U>3 = LO± 
^6 = ^7 = ^8 



(7) 



where the monomer weight and the corner weights are 
squared due to the fact that we are dealing with a pair 
of Majorana fermions glued together [l8| ■ 

Let us now turn to the description of the new method 
to efficiently simulate any vertex model in arbitrary di- 
mensions with generic (positive) weights a;, , including the 
fermionic models discussed above. For illustrative pur- 
pose we restrict the discussion to the 8-vertex model. 
The method is an extension of the so-called worm algo- 
rithm by Prokof 'ev and Svistunov Q . The configuration 
space of closed loops is enlarged to contain also open 
strings. For the GN model such an open string with ends 
at x and y corresponds to the insertion of two Majorana 
fields at positions x and y which is simply the Majorana 
fermion propagator 



G(x,y) 



Vie 



-Sgn 



(8) 



Similar interpretations of the open string can be obtained 
for other vertex models. The open string is now the ba- 
sis for a Monte Carlo algorithm which samples directly 
the space of 2-point correlation functions instead of the 
standard configuration space. This is the reason why the 
algorithm is capable of beating critical slowing down as 
demonstrated below: at a critical point where the corre- 
lation length grows large, the configurations are updated 
equally well on all length scales up to a scale of the order 
of the correlation length. 

In the vertex language the insertions correspond to the 
new vertex configurations depicted in the bottom row of 
Fig. [TJ A configuration containing a single open string 
corresponds to a loop configuration with two instances 
of vertex 9-16 which are connected by a string. Note 
that vertices 13-16, while present in the generic extended 
vertex-model, do not have a physical interpretation in 
terms of fermionic fields since they are explicitely for- 
bidden by Pauli's exclusion principle (fermionic lines are 
not allowed to intersect). Nevertheless they can also be 
included in the fermionic models, simply for algorithmic 
efficiency, and we do so in our implementation. 

The algorithm now proceeds by locally updating the 
ends of the open string using a simple Metropolis or heat 
bath step according to the weights of the corresponding 
2-point function. When one end is shifted from, say, x 
to one of its neighbouring points y, a dimer on the cor- 
responding bond is destroyed or created depending on 
whether the bond is occupied or not. In the process, the 




FIG. 2: Results for the ratios Zc i:j /Z with i,j = 0, 1 for the 
Majorana GN model on a 128 2 lattice as a function of the bare 
mass m. Dashed lines are the exact results calculated from 
the Pfaffians. Note that the curves and the data for Zc 10 /Z 
and Zc 01 /Z lie on top of each other. The inset shows the 
partition function ratio Z ( ^°/Z and illustrates how the zero 
mode at m = is reproduced. 



two vertices at x and y are changed from v 
and the move is accepted with probability 



x i u y 



to vL 



P(x — ► y) = min 



UJ„i 0J v i 



(9) 



in order to satisfy detailed balance. So a global update 
results from a sequence of local moves, and in this sense 
it is similar in spirit to the loop cluster update suggested 
in 0. 

Whenever the two ends of the open string meet, a 
new closed loop is formed and the new configuration con- 
tributes to the original partition function Z in one of the 
classes £oo, £iO)£oij £n- In this way the overall normal- 
isation is ensured, and expectation values can be calcu- 
lated as usual. From here it also becomes clear that the 
algorithm switches between the topological sectors with 
ease: as the string evolves it can grow or shrink in any 
direction and wrap around the torus as many times as it 
likes. Effectively, the algorithm simulates a system with 
fluctuating topological boundary conditions. 

In principle, the weight of the open string can be cho- 
sen arbitrarily, but the physical interpretation given by 
cq. (|HJ) suggests to choose the weights ujg to uie such that 
the open string configurations sample directly the 2-point 
correlation function. The open string configurations then 
provide an improved estimator for the 2-point correla- 
tion function [19j]. During the simulation one simply up- 
dates a table for G(x,y) as the string endpoints move 
around and the expectation value is obtained by form- 
ing (G(x,y))z = G(x,y)/Z. For the fermionic models 
we also need to keep track of the Dirac structure asso- 
ciated with G(x,y). This is most easily done by adding 
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FIG. 3: Results for the ratio Z® /Z in the Schwinger model at 
strong coupling for various system sizes. The inset shows the 
determination of the critical mass m c — —0.686506(27) while 
the main plot shows the collapse from finite size scaling con- 
sistent with a second order phase transition in the universality 
class of the Ising model (with critical exponent v = 1). 



FIG. 4: Integrated autocorrelation time ta of the energy for 
the Schwinger model in the strong coupling limit as a function 
of the system size L at the critical point m = m c . The line 
is a fit ta oc L z yielding z = 0.25(2). The inset shows the 
logarithmic dependence —18(3) + 7.5(6) ln(L) from fitting to 
L > 64. 



the product of the Dirac projectors along the string £, 
i.e. Il^ef P{^)i as a contribution at each step. Care has 
to be taken when the open string winds an odd times 
around a boundary on which we want to impose antiperi- 
odic boundary conditions for the fermions. In that case 
we need to account for an additional minus sign in the 
contribution to G(x,y). For the fcrmionic models where 
vertices 13-16 have no physical meaning, the weights W13 
to a; 16 can be tuned for algorithmic efficiency and do 
not follow any physically inspired rule. A good choice 
is to use the geometric mean of the weights ioi of those 
vertices that can be reached in one further update step, 
e.g. W13 = (UJ4UJ6UJ7) 1 / 3 . Finally, let us emphasise again 
that the algorithm described here is applicable to any 
vertex model, also in higher dimensions, as long as the 
weights are positive definite in well defined configuration 
classes. 

We have performed extensive tests of our algorithm 
by comparing to exact results known from Pfaffians (for 
the Major ana GN model) or from explicit calculations 
on small lattices. Simple observables are linear combina- 
tions of partition functions and ratios thereof, e.g. Zc i:j /Z 
with i, j = 0, 1. In Fig. [2] we show the results for the ra- 
tios Zc i3 /Z in the Majorana GN model on a 128 2 lattice 
as a function of the bare mass m. Dashed lines are the 
exact results calculated from the Pfaffians. Note that all 
partition function ratios are obtained in the same simu- 
lation. In the inset we also show the ratio Z®° /Z where 
Z®° = Zc aa — Zc al — Zcw ~ Zen is the partition function 
with fcrmionic b.c. periodic in space and time direction. 
In that situation the Majorana Dirac operator has a zero 
mode at m = (and at m = —2) and the system is 
critical. The inset in Fig. [5] illustrates that the algo- 



rithm can reproduce this zero mode without problems 
and that we can in fact simulate directly at the critical 
point. Conversely, we can use Z®°/Z = as a defini- 
tion of the critical point m = m c . In Fig. [3] we show 
our results for Z9 jZ as a function of the bare mass m 
in the Schwinger model in the strong coupling limit for 
various volumes. The critical point can be determined 
accurately with very little computational effort and we 
obtain m c — —0.686506(27) (cf. inset in Fig. [3]) from our 
simulations on the largest lattice with L = 512. Further 
improvement could be achieved by employing standard 
reweighting techniques as done in (l7| where they ob- 
tained m c = —0.6859(4). These calculations indicated 
a second order phase transition in the universality class 
of the Ising model (with critical exponent v ~ 1). Our 
results in Fig. [3] now confirm this by demonstrating that 
the partition function ratios Z9 jZ as a function of the 
rescaled mass (rn — m c )L v with v = 1 beautifully collapse 
onto a universal scaling curve. 

The efficiency of the algorithm and the fact that criti- 
cal slowing down is essentially absent is demonstrated in 
Fig. 2] where we show the integrated autocorrelation time 
ta of the energy as a function of the linear system size 
L at the critical point m = m c . (Similar plots can be 
obtained for the Majorana GN model.) The functional 
dependence on L can be well fitted (x 2 /dof = 1.28) by 
ta oc L z all the way down to our smallest system size 
L = 8. We obtain z = 0.25(2) which is consistent with 
just using the largest two system sizes. It is an amazing 
result that our local Metropolis-type update appears to 
have a dynamical critical exponent close to zero. The au- 
tocorrelation time may also depend logarithmically on L 
and a fit to L > 64 yields -13.8(1.9) + 6.6(4) ln(L) with 
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X 2 /dof = 1.00. 

In conclusion, we have presented a new type of algo- 
rithm for generic vertex models. It relies on sampling di- 
rectly 2-point correlation functions and essentially elimi- 
nates critical slowing down. We have successfully tested 
our algorithm on the Majorana GN model and on the 
Schwinger model in the strong coupling limit and found 
remarkably small dynamical critical exponents. The al- 
gorithm definitely opens the way to simulate efficiently 
generic vertex models (with positive weights) in arbitrary 
dimensions, in particular the GN model with any num- 
ber of flavours, the Thirring model, the Schwinger model 
in the strong coupling limit (in arbitrary dimensions) , as 
well as fermionic models with Yukawa-type scalar inter- 
actions, all with Wilson fermions. 
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